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We study small oscillations of the order parameter in weakly and strongly paired superconductors 
driven slightly out of equilibrium, in the collisionless approximation. While it was known for quite 
some time that the amplitude of the oscillations in a weakly paired superconductor decays as i -1 / 2 , 
we show that in a superconductor sufficiently strongly paired so that its fermions form bound states 
usually referred to as molecules, these oscillations decay as i~ 3//2 . The transition between these two 
regimes happens when the chemical potential of the superconductor vanishes, thus the behavior of 
the oscillations can be used to distinguish weakly and strongly paired superconductors. These results 
are obtained in the mean field approximation which may not be reliable in the crossover region 
between the strong and weak pairing, so we also obtain identical results within the two-channel 
model, which can be tuned to be reliable throughout the entire crossover, although it then describes 
a special type of interactions between the fermions which may be difficult to observe experimentally. 
Finally, we interpret the result in the strongly paired superconductor as the probability of the 
molecular decay as a function of time. 



PACS numbers: 05.30.Fk, 03.75.Kk, 03.75. Ss 

The study of quantum quenches, or the evolution of a 
quantum system when its parameters suddenly change, 
acquired prominence due to recent studies in Refs. OH]- 
Yet a particular case of it, the study of small oscillations 
of the order parameter in a perturbed superconductor, 
started more than 30 years ago in Ref. [3]. This paper 
showed that if a small perturbation is applied to a su- 
perconductor leading to a deviation of its gap from its 
equilibrium value, then this perturbation will evolve as 
cos(2Ao£) /Vt, where Ao is the equilibrium value of the 
gap. This work was continued in Ref. [4]. These authors 
showed that large oscillations of the order parameter in 
a superconductor do not decay but in fact continue for- 
ever (as long as the collisionless approximation remains 
valid) . Subsequently with the help of the exact integra- 
bility of this problem [5] it was shown that the exist a 
critical strength of the perturbation such that if a pertur- 
bation beyond that strength is applied, it excites oscil- 
lations which continue forever while perturbations below 
that strength still decay, Ref. [BJ [7] . 

At the same time, advances in atomic physics allow 
now to create superconductors out of ultracold fcrmionic 
atoms whose interactions can be externally controlled. If 
the interactions are weak, these artificial superconductors 
behave just like the usual superconductors of condensed 
matter physics (termed BCS superconductors). But the 
interactions can also be adjusted to be strong, in which 
case these superconductors become more akin to a Bose 
condensate (BEC) of diatomic bosonic molecules. As the 
interactions are tuned, the superconductor is said to un- 
dergo a BCS-BEC crossover (which is accompanied by 
the chemical potential changing from positive to negative 
values). The possibility of this crossover was discussed 
in a number of papers throughout the last four decades 
(see Refs. [8] EJ HQ]), and it finally was observed a few 
years ago (see Refs. [TTIIT2"] ). 



The unprecedented control over the interactions that 
the ultracold gases provide allows to change them quickly 
thus easily creating an initial out-of-equilibrium pertur- 
bation by a sudden change of their strength. A natural 
question which arises in this regard concerns the fate of 
the oscillations, decaying and persistent, as the supercon- 
ductor undergoes the BCS-BEC crossover. Refs. [131 E] 
explored the large nondecaying oscillations of a tunable 
superconductor using the ansatz of Ref. [3]. However, 
this method was found to break down as the supercon- 
ductor is tuned to the BEC regime, and it was not pos- 
sible to tell whether this was an artifact of the technique 
or an indication that an oscillating solution is indeed not 
possible in a strongly paired superconductor (see, how- 
ever, Ref. [T5]). 

To remedy this situation, in this paper we consider 
small oscillations of the order parameter in a tunable 
superconductor. We show that as long as the oscilla- 
tions remain small, they always decay regardless of the 
strength of interactions in the superconductor. However, 
while in the BCS superconductor the amplitude of these 
oscillations decay as l/Vi, in the BEC regime they de- 
cay as 1/i 3 / 2 . The transition between the two regimes 
happens exactly where the chemical potential fi is equal 
to zero. 

The theory we develop is based on the mean field 
(collisionless) approximation. Thus one can question its 
validity near the unitary point |16j which lies between 
the BCS and BEC regimes of the superconductors (at 
positive chemical potential) and is known to go beyond 
the applicability of the mean field theory [T7] (see also 
Ref. [TS] for recent studies of the dynamics near this 
point). Therefore, we repeat the calculation in the two 
channel model [19] describing BCS-BEC crossover with 
narrow Feshbach resonance [20] where mean field the- 
ory is applicable throughout and recover the same result. 
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This shows that an example of the BCS-BEC crossover 
exists whose order parameter oscillations unambiguously 
obey the scenario discussed here. 

In the "deep BEC" side of the crossover, the super- 
conductor becomes the Bose-Einstein condensate of the 
diatomic bosonic molecules. In this regime, the result 
l/£ 3 / 2 has a very simple interpretation. Indeed, as we 
will see in this paper, the order parameter of such su- 
perconductor is proportional to the probability of find- 
ing two fermions of opposite spins at the same position 
in space. As the superconductor is perturbed, the wave 
function of a pair of fermions no longer coincides with 
their bound state, but rather consists of the linear com- 
bination of a bound and excited states. The part of the 
wave function in the excited states moves off to infinity so 
that the probability amplitude to find two particles in the 
same spot decreases as l/t 3 ^ 2 reaching a limiting finite 
value at large t which in turn follows from the behavior 
of a three dimensional propagator of a free particle. In 
fact, in this deep BEC regime this is true not only with a 
small but also with a relatively large initial perturbation 
(the precise criteria are developed below), in which case 
this picture makes it clear that the order parameter de- 
cays to a value smaller than the equilibrium value, just 
like in the BCS regime with certain perturbations [7]. 
However, to fully analyze the case of large perturbations 
to the superconductor and to see whether nondecaying 
solutions are possible in the BEC regime, one needs to 
take advantage of the integrability of the equations of 
motion, something which was done successfully for the 
weakly paired superconductor and whose strongly paired 
superconductor applications are left for future work. 

In what follows, we describe the derivation of these 
results. Consider spin- 1/2 fermions interacting via a 
short range attractive s-wave interaction of strength A. 
Within the mean field approximation, these fermions 
can be studied on a strictly classical level. We ac- 
complish this by introducing the Cooper pair number 

n p = ^aj p aTp + ^\p a lp X ) T ne evolution of this num- 
ber obeys the following classical Hamiltonian, a classical 
version of the Anderson- Richardson (or reduced BCS) 
Hamiltonian describing these kind of fermions, 



H = 2e p n p - - ^ 



p.q 



^Jn p (l - n p )n q (l - n q ) cos(4> p - 4> q ), 

(1) 

where e p — p 2 /(2m) is the free fermions dispersion, 4> p 
is the phase variable canonically conjugate to n p , and 
V is the space volume. The equations of motion of this 
Hamiltonian have the following solution, corresponding 
to a stationary superconductor 



-2(it, n° p 



1 



1 
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where as always /i is the chemical potential, t is time, 



£p = 6 p~^j E p — y £, p + Ag, and finally Ao is the equilib- 
rium gap of the superconductor which can be found from 
the gap equation 1 = ^7 ^2 p -jjr- ■ Now consider the initial 
conditions for the motion described by Eq. ([IJ consisting 
of a small perturbation to the stationary solution Eq. Q . 
In fact, it is most physical to take as initial conditions 
Eq. ([2| for a slightly different value of A = A + <5A 
as well as the slightly different value of fi (corresponding 
to a superconductor whose interactions were slightly per- 
turbed in the initial moment of time) . Then it is straight- 
forward to see, with the help of the particle conservation 
condition ^ p 5n p = 0, that the initial conditions read 



ip = n° + 5n° p , dn° = I £ p - 



fo 

h 
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where /o and fo are defined below in Eq. (|9|. We then 
expand the Hamiltonian Eq. (|T|) about the stationary so- 
lution to obtain the quadratic Hamiltonian for the devia- 
tions 8n p , S(f>p. Subsequently we construct the solution to 
the equations of motion of this Hamiltonian by using the 
method of Green's functions. The expanded Hamiltonian 
takes the form 



SH = 



2^1 



f)p*§* pqt 



b q + 5n p K pq Sn q ] , 
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This represents a collection of harmonic oscillators, la- 
belled by the index p. Now it is possible to construct a 
retarded Green's function corresponding to these oscilla- 
tors, with the end result 



Sn p (t) =i^2 J 
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where 

Gpq(fi) 



E 2 
p 
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a / (n%g q + e 2 (n 2 ~a 2 ))- f (t P E 2 + £ q n 2 ) 
e 2 ) (n 2 - ei) \f 2 + f 2 (A 2 - n 2 ) 
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Here we introduced the functions 

A 1 ~ A 



/ = 



2V 



E 



E p (tt 2 
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In what follows we will also need these same functions 
evaluated at f2 = 0, or 



fo = - 



A 

2~V 
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The derivation of Eq. ^ is technical and not particularly 
instructive, so it will be published elsewhere. 

Armed by the explicit expression for the fluctuations 
Sn p (t) we can now discuss how to calculate the time de- 
pendent gap function. It is easiest to study the square of 
its absolute value, 



A 2 

V 2 



^2 V U P( 1 ~ n p) n q( l - U q) COS(0 P - (j> q ). 

(10) 

In turn, this quantity can be decomposed into the sum 
of A 2 ,, the square of the unperturbed gap, and the per- 
turbation S |A| 2 , given by 



2 A 



5\A(t)\ 2 = v Y,t P 5n p (t). 



(11) 



We combine this with the initial conditions Eq. Q as 
well as with Eqs. ^ and Q to find that 



S\A(t)\ 



2 iA SA Q 
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ndn e - 2mt x (12) 



' — 00+2O 

~P + P (A 2 -fi 2 )-// A 2 -// 2 // 



O 2 



p + p (A 2 - CI 2 ) 



This equation represents the main result of this paper. 
White reducing to the main result of Ref. [3] at weak 
pairing, it represents the generalization of their result to 
a superconductor of an arbitrary pairing strength. The 
integration here goes over the straight line slightly above 
the real axis, as shown on Fig. [T|with a dashed line. We 

1 Q 



FIG. 1: The contour of integration in Eq. ( 12 1 in the complex 



plane of Cl. Dashed line shows the initial contour, while solid 
black line, together with a dashed-dotted semicircle, is the 
deformed contour. The dashed-dotted part of the new contour 
can be deformed away to infinity and does not contribute. 

expect the functions / and / to have singularities on the 
real axis where |f2| hits the minimum of E p and cuts at 
larger |f2|. It may seem that there is also a singularity 
at CI — 0, but it is actually fictitious as the numerator 

At t < the contour 



of Eq. (12 1 vanishes at CI 



0. 



can be closed in the upper half plane resulting in zero 
for the integral as expected. At positive time t > we 
can deform the contour as shown on Fig [l] by a solid 
line with arrows. Then at large positive times, the main 
contribution to the integral comes from the vicinity of the 
singular points, or the turning points of the new contour. 



To evaluate / and /, we replace summation over p 
in Eq. Q by the integration over Vd 3 p/(2ir) 3 . In the 
BCS regime, it is standard to pass to the integration over 
the variable £ = e p — p so that d 3 p/(2n) 3 ss vd^, where 
v = \Jvn?iil(j{ 1 \pi) is the density of states, and then 
extend the integration over £ all the way to the entire 
real axis. Then / = as its integral is antisymmetric in 
£. At the same time, / is given by 



f = i 



ci^/a 2 - ci 2 



/ ci 2 ci 



(13) 



It is then easy to see that the integrand in Eq. ( 12 1 has a 



singularity as Cl approaches A which goes as 1 /y/Cl — A 
in agreement with the discussion above (Ao coincides 
with the minimum of E p ). Introducing the variable 
s = Cl — Aq and taking into account that at t 1/Ao 
only small s contributes to the integral, we find after 
adding both left and right contours (which are complex 
conjugate of each other) 



*iA(*)r 
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Alternatively, in the BEC regime where [L < 0, the inte- 
gration measure in Eq. (JsJ) must be kept as Vd 3 p/(2ir) 3 
and both functions / and / are not zero. The minimum 
of E p is equal to i^min = \J + Aq. The functions / 
and / are now finite when |fi| reaches -E m i n , but they 
are not regular at that point. Generally they behave as 
/ ~ const + i^/il — E min . The origin of such a drastically 
different behavior lies in the fact that the minimum of 
E p is reached when p = 0. Then, the potential diver- 
gence in the integral at = E m i n is removed by the p 2 
coming from the measure of integration p 2 dp. All this 
can be checked explicitly in the deep BEC regime where 
3> A and the analytic expressions behaving precisely 
in this way can be obtained by direct integration. 

Now if we substitute this into Eq. ( [12] ) we find that the 
constant piece is single valued as CI goes from the upper 
to lower branch of the contour on Fig. [T| and therefore 
its contribution from each of the contours cancels. So it 
is the square root part which contributes. This gives at 
large times t > 1 /E m i n 
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(15) 



As we see, the transition from the behavior Eq. (14) to 



Eq. ( 15 1 occurs exactly at /j, — as this is the point where 
the minimum of the spectrum E p shifts to p = 0. 

One can worry if this derivation really works for a su- 
perconductor in the vicinity of the unitary point where 
the mean field theory breaks down. To have a control- 
lable theory, we repeated this calculation for the two 
channel model, based on the equations of motion intro- 
duced in Ref. |13j . That model has an additional coupling 
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g such that if g 2 m 2 /n 1 / 3 ^> 1 (here n is the density), 
then the two channel model is equivalent to Eq. (JT|) con- 
sidered earlier in this paper. In the opposite limit, the 
two channel model still undergoes a crossover but with 
the mean field theory applicable throughout. The answer 
for the two channel model can be worked out using the 
same methods as the ones described here. It looks almost 



identical to Eq. ( 12 1, except in all occurences of / and /o 
one needs to substitute / — > / — 2 A/5 2 , f — > fo — 2 A/5 2 . 
This replacement does not change any of the arguments 
presented here for the one channel model, thus all the 
conclusions remain valid. 

Finally, we observe that in the deep BEC regime all 
n° <C 1. Then the Hamiltonian Eq. (Jl| can be signifi- 
cantly simplified, reducing in this limit to the Schrodingcr 
equation of one pair of fermions in a delta-function po- 
tential, as discussed at length in Ref. [ID]. Replacing 
1 — rtp pa 1, and introducing the pair wave function 



— V"'' 
ing to 



e^*>, we find H = 2e p \tp p 



Dp A 



id t ip = 2e p i> p 



1 



lead- 



(16) 



The "gap function" in this language is A = A J2 P *Pp/V 
which has an obvious meaning of the probability of find- 
ing two opposite-spin fermions at the same point in space. 



We can solve the linear equation (16) with arbitrary ini 



tial conditions directly. The solution for the gap function 
reads 



A(t) 



n = 
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V 
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with ipp being the initial value of ip p . To compute the 
integral over f2, we deform the contour of integration in 
a way similar to that shown on Fig. [I) Unlike Eq. (jl2|, 
the integrand here has a simple pole at a negative value 
of Cl = fi& where H(£lb) — — 1, but similarly to Eq. (12) 
in the BEC regime, it has cut at £1 > where it behaves 
as const + iy/Q. The pole corresponds to the bound state 
and results in a contribution to A(t) corresponding to 
the bound molecule (that is, the equilibrium solution). 
The cut produces the decaying behavior going precisely 
as 1/t 3 / 2 , thus overall |A(i)| 2 « const+cos(fi b £)/t 3 / 2 (for 
comparison with Eq. ( 15 1 observe that fib = 2\x in this 



regime). However now we see that this decay is correct 
for arbitrary initial conditions, as long as the applicabil- 
ity conditions \i/j p \ -C 1 are satisfied. The value |A(t)| 2 
decays to is smaller than the equilibrium value of the gap 
function, since it is reduced by the overlap, smaller than 
1, of the initial state and the bound state. The only ini- 
tial distributions out of reach of this method are the ones 
where | ip p | are no longer much smaller than 1 . 



In conclusion, we showed that the amplitude of the 
small oscillations of the order parameter in a supercon- 
ductor decays as l/i 1/2 in the BCS (/z > 0) and as l/i 3/2 
in the BEC (/i < 0) regimes. Applicability of this result 
to the Fermi gas close to the unitary point remains an 
open question. It would also be very interesting to see if, 
with modern experimental techniques, these oscillations 
can be detected in real superconductors and whether this 
theory can be generalized to superconductors with the d- 
wave symmetry of the gap. 
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